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Molecular Dynamics simulations of water confined in nanometer sized, hydrophobic 
channels show that water forms localized cavities for pore diameter > 2.0 nm. The 
cavities present non-spherical shape and lay preferentially adjacent to the confining 
wall inducing a peculiar form to the liquid exposed surface. The regime of localized 
cavitation appears to be correlated with the formation of a vapor layer, as predicted 
by the Lum- Chandler- Weeks theory, implying partial filling of the pore. 

I. INTRODUCTION 

Nanotubes [1] are important building blocks of nanocomposite materials and nanomachin- 
ery. When immersed in ionic solutions, nanometer-sized pores can be used for the detection 
and analysis of electrolytes and charged polymers, such as DNA [2], and to understand 
ion transport and selectivity in biological channels A full understanding of the phase 
behavior of confined water is a pre-requisite to interpret adsorption and conductivity data. 
In particular, experiments and simulations have shown that, for channels of sub-nanometer 
radius water exhibits an intermittent filling/drying transition resulting from an underlying 
bi-stable free-energy landscape separated by a thermally activated barrier, a behavior spe- 
cific to water and not observed in a monoatomic fiuid [HIS]. The switch between empty /filled 
states upon ion translocation has been advocated as a gating mechanism in passive biological 
channels [HI E]. 
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Cavitation of water under confinement is related to the long-ranged attractive forces ex- 
erted by hydrophobic bodies. Two alternative mechanisms have been proposed for explaining 
the occurrence of these forces: on one hand, the presence of air-filled nanobubbles at the 
solid surface drives the attraction via a bubble-driven mechanism; on the other hand, long- 
range correlations in the critical region or near the spinodal line could induce an extended 
density depletion at the surfaces which results in the attraction P [91 [TO]. 

The first fine of reasoning relies on the pre-existence of air-filled cavities at the solid surface 
[TT] . However, the mechanism for their stabihzation is not fully understood since the Laplace- 
Young equation prescribes an internal pressure of the order of 10^ atm for nanometer-sized 
bubbles [12]. One question pertains how such large internal pressure can sustain stable 
or metastable cavities, and if these are independently nucleated at the sohd, corrugated 
surface. It has been suggested that line tension can act to stabilize such small cavities |T3] . 
by reducing the curvature of the bubble base on the solid surface. As a matter of fact, 
micron-sized air bubbles are commonly observed in proximity of hydrophobic corrugated 
surfaces and, in the recent literature, there is growing evidence for the presence of bubbles 
on corrugated surfaces at the nanoscale. In particular, recent AFM experiments on silicon 
oxyde wafer surfaces of controlled roughness reported on bubble formation at the solid-water 
interface [H], while syncrotron x-ray reflectivity measurements reported on a small depletion 
layer in lieu of localized cavitation |T6] . 

The second interpretation is based on the formation of a vapor layer in contact with the 
solid surface. Such depletion layer would involve a high entropic cost, growing with the 
extension of the exposed surface, but compensated by volume-dependent forces arising from 
molecular reorganization. The Lum-Chandler- Weeks (LCW) theory formulates the hy- 
drophobic effect in microscopic terms based on the competition between interfacial and bulk 
forces which therefore depend on the solute surface/volume ratio. The molecular mechanism 
underlying the force balance is commonly ascribed to the distortion of the hydrogen bond 
network. In particular, for small spherical solutes water arranges in a clathrate structure, 
whfle for larger size the ensemble of hydrogen bond vectors, 0-H - H, points preferentially 
towards the solute. Therefore, if the curvature of the embedded body is low, the distortion 
is large and entropic effects prevail over the enthalpic ones, and vice versa for small solutes. 

The cylindrical pore represents an opposite case to the spherical solute, where now the 
conflning surface is concave. However, a similar enthalpic/entropic competition is expected 
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and should depend on the pore diameter. Thus, a relevant information pertains if, for 
pore diameter larger than that characteristic of intermittent behavior, the liquid phase 
is completely stabihzed or if some peculiar behavior still takes place. The occurrence of 
cavitation in water confined by cylindrical pores has been recently discussed in view of 
understanding the hysteresis involved in capillary evaporation [T8] . 

We report here the results of Molecular Dynamics (MD) simulations undertaken on cylin- 
drical, hydrophobic pores of finite length. In section II the simulation set-up and numerical 
details are described and in Section III the Molecular Dynamics data are illustrated. In 
Section IV the LCW theory is described and the numerical solutions are compared with the 
simulation results. The last Section draws some conclusions and analyzes recent experimen- 
tal observations. 

II. SIMULATION DETAILS 

The geometry of the simulation set-up is depicted in Fig. [T] The two pore mouths are 
connected to the same, periodically folded, reservoir of quasi-bulk water. The pore diameter 
d and pore length L are varied to the following pairs of values {d, L) = (1.5,4.0), (2.0,4.0), 
(2.5,4.0), (2.5,8.0), (3.0,4.0) and (3.5, 4.0) nm, as used to label the runs. Simulations be- 
come rapidly costly with the pore size and we did not attempt to extend our simulations 
to larger geometries. Besides the nominal values, the effective cylinder diameter and length 
were evaluated via the criterion that the repulsive wall-oxygen potential be less than ksT, 
giving a reduction of the effective diameter and increase of effective length by 0.46 nm [5]. 

The simulation box is periodic in all three directions with dimensions depending on the 
pore geometry. The number of water molecules for all simulations ranges between 852 and 
2464 units. Water is represented via the SPC/E computational model [22]- The confining 
oxygen-wall potential is modeled as a smooth surface generated by carbon atoms distributed 
uniformly over the surface. The wall-water interaction acts between each water oxygen and 
a smooth Lennard-Jones potential integrated over the dark region of the wall in Fig. [TJ 
Each carbon atom carries a Lennard-Jones potential with parameters a = 0.345 nm and 
e = 0.7294 A; Jmo/"^ p] and standard Lorentz-Berthelot rules are used to construct oxygen- 
wall interactions. Electrostatic interactions between charged oxygen and hydrogen atoms 
are computed with the Ewald method via the Smooth Particle Mesh Ewald implementation 
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|19| . Moreover, the confining medium is taken as non polarizable. 

Special care is taken to keep the reservoir density under control, such that upon empty- 
ing/filling of the channel the reservoir maintains a bulk-like behavior with an average mass 
density of Igr/cm^ away from the pore. This is obtained by varying over time the length 
of the reservoir, in direction parallel to the pore axis, by a Berendsen type of piston |20] . 
with a characteristic couphng time of 10 ps, during the equilibration and production runs. 
The feedback sets the average density in a stripe of thickness t = 10 nm placed at distance 
/ = 7nm away from the pore (see Fig. [T]). The system temperature is controlled via a Nose- 
Hoover thermostat [21] in order to avoid anomalous temperature drifts during both the 
equihbration process and the subsequent production runs. The system is simulated at 300 k 
while some auxiliary simulations are made at 280 and 320 K. The evolution of the system is 
followed over times of 100 (equilibration) and, subsequently, of 1 ns (production). Within 
the simulation time window, stationarity is monitored by following the number of molecules 
populating the channel and the formation and size distribution of cavities. In order to moni- 
tor the effective stationarity of the equihbrated system, the run of the (2.5, 4.0) nm geometry 
has been further extended to 3 ns, without observing any departure from stationarity. 

In order to detect the cavities we have used a coarse-graining procedure by tesselating 
the space with cubic cells of edge 0.02 nm. An empty cell is defined by having the distance 
from any oxygen atom greater than 0.3 nm. The wall position is defined by the largest radial 
distance of an oxygen atom from the pore axis and adding an offset of 0.1 nm. A cavity is 
defined as the cluster of continguous empty cells which do not belong to the wall region. 
The distinction among clusters is made via a graph algorithm and sorted according to size. 
In this way, the identity of cavities is clearly established. 

III. MOLECULAR DYNAMICS RESULTS 

The simulation data present two distinct behaviors: for the geometries (1.5, 4.0) nm at 
T = 300 if and (2.5, 4.0) nm at T = 280 we observe emptying of the pore, with a large 
vapor region extending between the two pore mouths; on the other hand, for the geometries 
(2.0,4.0), (2.5,4.0), (2.5,8.0), (3.0,4.0) and (3.5, 4.0) nm, water persists inside the pore at 
density close to liquid state, with formation of nanobubbles in proximity of the confining 
wall. 
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At first we report on the process of emptying of the channel, as depicted in Fig. [2} Given 
an initial, apparently stable hquid inside the channel, a number of distinct lateral cavities 
are formed in proximity of the cylinder wall. Subsequently, one of such bubbles grows in 
size, rapidly extending over the whole pore volume. The final state appears to be stable 
over the 1 ns time scale. The homogeneous growth of the vapor region is consistent with 
the recently proposed model of capillary evaporation in hydrophobic pores [T8] . 

Intermittent emptying/filling of the channel has been previously reported at T = 300 K 
for d < 1.5 nm and L < 1.0 nm (and eventually for larger diameters at lower temperatures) 
by using MD simulations with the same computational model as used here [5]. In the 
previous work, it was found that the time scale of empty ing/filhng oscillations is of the 
order of ~ 100 ps. In the current case, the direct observation of intermittency is prohibitive 
in terms of CPU time. In fact, given the larger pore diameters, intermittency could take 
place over timescales longer by one order of magnitude, or more [7]. Therefore, we preferred 
to focus on the interfacial structuring of the liquid in the filled state. 

It is interesting to note that the geometry (2.4, 4.0) nm exhibits pore emptying for T = 
280 K and pore filling for T = 300 K. The fact that emptying/filling of the pore is sensitive 
to differences in temperature as small as AT = 20 is an indication that the underlying 
bulk phase diagram, with the rather narrow region of bulk liquid water and the close-by 
hquid-vapor coexistence, is an important driving factor for the confined fluid. 

For the filled channel, the liquid interface displays the patterns shown in Figs. [3] and 
[4] The interface is rather corrugated by the presence of cavities appearing in proxim- 
ity of the confining wall (lateral cavities), and rarely observed as spherical shape located 
around the cyhnder axis. The cavities resemble spherical caps or rounded cones, with 
volumes much larger than the molecular size {vw — 0.027 nm^). The average volumes 
range between = 5.7 for the {d, L) = (2.0, 4.0) nm geometry and {v)/v.uj = 21.0 

for {d,L) = (3.5, 4.0) nm. The peculiar pattern of the liquid exposed surface is rather rep- 
resentative of the highly cohesive character of water. It should be noticed that the shape 
of the vapor cavities along the transversal section does not appear to be flattened against 
the wall, while in the longitudinal direction visual analysis does not allow to draw a clear 
conclusion. We infer that line tension effects can be important only along the longitudinal 
direction but not on the transverse direction. In terms of dynamics, the cavities present 
high mobility with fast reshaping on the Ips time scale. 
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The distribution of cavity size is illustrated in the normalized histograms of Fig. [5| A 
systematic increase of cavities of larger volume with the pore diameter is visible, accompanied 
by a depletion of smaller, molecular-sized ones. The distribution becomes more shifted to 
the right as the diameter of the cylinder grows. In contrast, water in proximity of a planar 
hydrophobic surface exhibits cavities with reduced volumes, basically indistinguishable from 
fluctuations of size of the molecular volume Vy^- 

In fig. [6] we report the same histograms but we filter out the contributions arising from 
cavities of volume v < v^^. By rescahng the abscissa by the cyhnder surface area, the his- 
tograms become peaked around the same value v/{vyjndL) ~ 2nm^^, and width increasing 
with the pore diameter. However, the integrated cavity volume per unit surface is mostly 
contributed by the right tafls of the distributions, which appear to be equal for the (2.5, 4.0) 
and (3.5, 4.0) geometries but significantly shifted to the left for the (1.5, 4.0) case. Therefore, 
cavitation does not seem to grow further as the diameter is larger than 2.5 nm. 

The longer channel does not appear to affect the histogram, indicating that the bubble 
formation does not depend on interfacial effects arising from the finite pore length. In 
principle, the finite length could affect the liquid/vapor balance and the drying transition. In 
fact, by using an elementary macroscopic argument [I2], the difl'erence in free-energy density 
between the liquid and vapor states is approximated by AQ/iidL = -^i^ — ^lydjlL where 
7;^ and 7^ are the liquid-wall and liquid-vapor surface tensions under ambient conditions, 
respectively, and bulk contributions to the free energy have been neglected. In other words, 
if finite channel effects are in place, by making the channel longer, one should move away 
from the vapor branch. By taking 7/^„ = 7 AksT /nrn? and 7^ = 2QkBT /nrn? it is seen that 
the geometry {d,L) = (2.5, 4.0) nm is close to the region where the free-energy difference 
changes sign [5]. Our simulation at room temperature for {d,L) = (2.5, 8.0) nm showed 
that the channel remains filled as much as the (2.5, 4.0) nm geometry, without appreciable 
differences in the distribution of cavity volumes. Conversely, the sensitivity to the pore length 
appears in the fiuctuations of the number of adsorbed water molecules, which approximately 
drops by a factor two in the longer channel. We interpret the weak size dependence of 
cavitation as due to the corrugated, largely exposed surface of the liquid with respect to the 
simple macroscopic model, that renders finite size effects negligible already at L = 4nm. 
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IV. LUM-CHANDLER- WEEKS PREDICTIONS 

In this section we compare our data to the predictions of the LCW theory for an infinitely 
long pore. Schematically, the theory can described as follows. Let us first consider the 
decomposition of the microscopic density n{r) into a slowly varying component, ns{r) , 
and a fast component, n(r) — ns(r). The LCW theory builds upon the well-known square- 
gradient theory for liquid-vapor coexistence [25] by taking into account self-consistently the 
fast oscillations in density. According to the square gradient approximation, the local grand 
potential —W{n) is related to the laplacian of the density via 
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mV'^nir) (1) 



where m is an effective parameter derived from the underlying interatomic potential f(r), 
-| drr*v{r) (Random Phase Approximation [25]). By applying a coarse graining 



procedure, the local density is replaced by the function n — n = n + -g-V n where A is 
the characteristic length of the coarse graining procedure. Therefore, the starting equation 
of the LCW theory is derived from eq.Q rewritten in terms of n, and asserting that this 
equations holds for the slow component n^, thus 

dW\ 2m, , , , 

for which we choose A = 0.3 nm and m = 919.9 x 10^^ Jnm^mo/"^. To implement eq. 
Q one needs prior knowledge of the bulk equation of state of water near the liquid-vapor 
coexistence. Following previous studies [23], we used a quartic dependence plus a linear 
correction term which models the proximity to coexistence, 

W{n) = a{n — riif [n — Ugf + I3{n — rii) (3) 

where a = 2.89 x 10"^^ Jnm^, (3 = 3.04 x 10~^^J and the liquid and vapor densities are 
rii = 32.94 nm"^ and Ug = 7.7 x 10""^ nm"^, respectively [23]. By solving the equations in 
cylindrical coordinates and forbidding any angular symmetry breaking, we did not attempt 
to observe the formation of cavities of given shape. This choice was motivated by the narrow 
range of numerical stability found during the solution of the LCW equations. 

The link between the slowly varying and the complete density fields is estabhshed within 
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the gaussian density fluctuations approximation, by writing [23] 

n(r) = ns{r) — J (ir'c(r')x(r, r') (4) 

where the response function x(r, r') = {6n{r)6n(r')) is approximated by r') ~ ns{r)6{r — 
r') + ns{r)ns{r')h{\r — r'|) and h{\r — r'|) — 1 is the bulk pair correlation function. Moreover, 
c(r) is the water- waU direct correlation function. Given ^^(r), eq.Q is solved to obtain c(r) 
and n(r) via the Nystrom numerical procedure [26]. The two coupled integro-differential 
equations ([2]) and are solved iteratively. We did not include the correction to the LCW 
theory due to the attractive part of the wall-oxygen potential [23]. Vice versa, the effective 
diameter of MD and the diameter used in LCW calculations, d, were considered as equivalent 
control parameters. 

We have found that for diameter d < 2.6 nm the LCW equations predict complete 
pore emptying, in agreement with previous simulation and experimental results, but be- 
low our MD data which exhibit emptying (or eventually intermittency) for effective diameter 
< 2.0 nm. However, in this geometry the LCW predictions do not distinguish between inter- 
mittent and cavitating behavior. In the range 2.6 < d < 3.8 nm the theory predicts partial 
fining of the pore, as fllustrated by the density profiles in Fig. 4, where the micro-phase 
coexistence is accompanied by a significant density depletion near the wall. Moreover, the 
LCW equations predict a gradual increase of the vapor layer as the diameter lowers, whilst 
the MD results show a maximum which varies slowly with the pore diameter. Following 
previous authors [23], we attribute this difference to the presence of the weak attractive 
tail in the water-wall potential. Finally, at larger diameters, the LCW curves show weakly 
structured profiles, with a characteristic non-wetting shape near contact. In essence, the 
LCW theory predicts a cross-over between empty and partially fiUed states for d = 2.6 nm, 
larger than the MD results where the cross-over appears in the 1.5 : 2.0 interval. This 
discrepancy may be attributed to the one-dimensional solution of the LCW equations and 
to the missing attractive tafl in our theoretical treatment, causes which artificially stabihze 
the vapor phase, or to the undetermined location of the confining surface, whose diameter 
can vary by about 0.5 nm. 

Notwithstanding the neghgible correlations emerging from the LCW solutions, we can 
compare MD and LCW data regarding the number of water molecules filling a channel of 
length 4 nm as a function of the channel effective diameter d. From Fig. [s] it is apparent 
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that the LCW model shows systematicaUy lower values than the MD data in the region 
where the theory exhibits partial fiUing. The MD data are intermediate between the LCW 
values and the simple model of uniform fiUing of the pore at the water liquid density. As 
for the discrepancy in the cross-over diameter, the lack of longitudinal symmetry breaking 
and the weak attractive tail in the simulated water-wall potential, might explain the smaller 
number of water molecules predicted by the LCW theory. 

V. CONCLUSIONS 

The behavior of water in contact with an extended hydrophobic surface has attracted 
considerable attention on account of the intrinsic thermodynamic problem, which has signif- 
icant imphcations for protein folding |23j|, for ionic transport in biological channels |E] and 
for the boundary conditions for fluid flow in microchannels [28] . 

Our study clearly demonstrates the spontaneous formation of vapor microcavities on 
the nanometer scale for water in contact with a concave cylindrical surface of diameter 
d > 2.0 nm. The cavities are locahzed both in the angular and, more importantly, in the 
longitudinal directions. To our knowledge, this surprising result is the first observation for 
water confined in such geometry. 

The results are particularly interesting for the on-going debate on the hydrophobic ef- 
fect, in which a number of different interpretations have been put forward, such as entropic 
effects due to molecular rearrangement, electrostatic effects and spontaneous cavitation due 
to the metastability of the fiuid jlO]. Recent experiments have focused on analyzing water in 
contact with a hydrophobic surface and found contrasting results. In particular, either cavi- 
tation [m [T3] or an extended low-density depletion layer [IH] have been observed. The fact 
that water around a concave, but smooth, hydrophobic surface forms short-ranged islands 
of vapor, modulated by the surface curvature, suggests that the contrasting experimental 
observations arise from the sensitive structural response of water to the roughness of the 
solid surface. 

Moreover, our study sheds some light on the experimental observation of hysteresis in 
water intrusion/extrusion cycles in pores and the interpretation based on homogenous nu- 
cleation [18]. The simulation data underhne the metastable character of highly confined 
water. However, the rich structural patterns exhibited by the interface are absent in a sim- 
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pie fluid and, therefore, are unlikely to be explained in terms of a macroscopic approach. 
When comparing the simulation data with a quasi-microscopic treatment, namely the Lum- 
Chandler- Weeks theory, the density profiles and the cross-over diameter between empty and 
partially filled states of the latter were found to agree only qualitatively with the simulation 
data. The discrepancies were explained on the basis of the one-dimensional solution of the 
LCW equations. 

Finally, we wish to comment on some recent measurements on ionic conductance in 
nanopores of diameter c? ~ 10 nm showing that transport displays an anomalous response, 
with a five orders of magnitude reduction in the current spectral density power and a strongly 
noisy response with respect to bulk behavior j29]. This has been attributed to the presence of 
spherically shaped air bubbles trapped inside the nanopores such that the translocating ions 
encounter a two-phase filled region. Although the differences in pore diameter between the 
presently simulated and the experimental systems might seem large, one might expect that 
water cavitation is still effective in the wider pores or that cavitation is locally enhanced 
by the translocating ion. The effect of such cavitation would induce non trivial wall-ion 
dielectric interactions, modulated by the imperfect screening of water, and non trivial hy- 
drodynamic forces. Such scenario would imply a coupled ion-bubble transport mediated by 
the microscopic hquid-vapor coexistence. 
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Figure Captions 

Figure 1: Schematic representation of a cut through the simulation cell. The dark area is the 
confining cylindrical wall and the shaded area is the region where water is maintained at mass 
density of Xgrjcvr? (see text for details). 

Figure 2: Snapshots of the (1.5, 4.0) nm channel at four different instants separated by lOps showing 
the process of pore emptying. The pictures are produced with the surface generation algorithm of 
the VMD software [21] by taking into account the oxygen atoms only. 

Figure 3: Snapshot of an instantaneous configuration of the (2.5, 4.0) nm channel, generated as for 
FiglU 

Figure 4: Evolution of cavities (shaded regions) formed in the (2.5, 4.0) nm pore at 300 and 
analyzed on a slab of thickness 0.1 nm placed at the midpoint of the pore axis. The four snapshots 
(a,b,c,d) correspond to successive configurations separated by \ps. The white region corresponds 
to the wall region, grey to water filled regions, black to the empty regions (see text for details). 

Figure 5: Frequency of occurrence of cavities with volume v normalized by the volume per water 
molecule {v^ = 0.027 nm^) for (d, L) equal to (2.0, 4.0) nm (solid line), (2.5, 4.0) nm (dashed line), 
(2.5, 8.0) nm (long dashed line), (3.5, 4.0) nm (dot-dashed line) and for water close to an infinite 
plane (dotted line). All data are at T = 300 ET. To highlight the dependence on the pore geometry, 
the normalized frequency of occurrence P{v) is multiplied by the volume. The inset displays the 
unsealed frequency (same symbols of the main plot). 

Figure 6: Histograms as in Fig. |5] but for cavity volumes larger than ■ The three curves refer 
to (d,-L) = (2.0,4.0), (2.5,4.0) and (3.5, 4.0) nm, with line styles as in Fig. |5] Abscissa have been 
divided by the cylinder surface area. 
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Figure 7: Radial density profiles divided by the density of liquid water n;. The radial coordinate 
is rescaled with the effective radius (MD data) and the value d/2 (LCW data). Upper panel: MD 
profiles at 300 A; for nominal pore diameter d = 2.0 nm (black), 2.5 nm (red), 3.0 nm (green) and 
3.5 nm (blue) and length L = Anm. Middle panel: LCW profiles for d equal to 2.8 nm (black), 
3.0 nm (red), 3.4 nm (green), 3.6 nm (blue). Lowner panel: LCW profiles for d equal to 3.8 nm 
(black), 4.0 nm (red), 6.0 nm (green), 10.0 nm (blue). 

Figure 8: Number of water molecules occupying channels of length 4nm obtained from MD (filled 
circles), LCW solutions (open squares) and a uniform filling of the channel (dashes curve). The 
vertical jump in the LCW curve signals the cross-over between empty and partially filled states. 
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